< 



Implicit-explicit Runge-Kutta schemes and applications to 
hyperbolic systems with relaxation* 

Lorenzo Pareschi^ Giovanni Russo* 

May 6, 2004 

O ' 

D . Abstract 

We consider new implicit-explicit (IMEX) Runge-Kutta methods for hyperbolic systems of 
conservation laws with stiff relaxation terms. The explicit part is treated by a strong-stability- 
preserving (SSP) scheme, and the implicit part is treated by an L-stable diagonally implicit 
Runge-Kutta methods (DIRK). The schemes proposed are asymptotic preserving (AP) in 
the zero relaxation limit. High accuracy in space is obtained by Weighted Essentially Non 
Oscillatory (WENO) reconstruction. After a description of the mathematical properties of 
f l ■ the schemes, several applications will be presented. 
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1 Introduction 

, Several physical phenomena of great importance for applications are described by stiff systems of 

C ' differential equations in the form 

(N . 

OV d t U = F{U) + -R(U), (1) 

where U = U(t) G R N , T, R : R N -> R N and e > is the stiffness parameter. 

System (JTJ may represent a system of N ODE's or a discretization of a system of PDE's, 
such as, for example, convection-diffusion equations, reaction-diffusion equations and hyperbolic 
■ systems with relaxation. 

In this work we consider the latter case, which in recent years has been a very active field of 
research, due to its great impact on applied sciences [T31|3T]. For example, we mention that hyper- 



bolic systems with relaxation appears in kinetic theory of rarefied gases, hydrodynamical models for 
semiconductors, viscoelasticity, multiphase flows and phase transitions, radiation hydrodynamics, 
traffic flows, shallow waters, etc. 

In one space dimension these systems have the form 

dtU + d x F(U) = -R(U), x G M, (2) 

£ 

that corresponds to (JXJ) for J-(U) = —d x F(U). In ([2j the Jacobian matrix F'(U) has real eigen- 
values and admits a basis of eigenvectors VUG M. N and e is called relaxation parameter. 

The development of efficient numerical schemes for such systems is challenging, since in many 
applications the relaxation time varies from values of order one to very small values if compared 
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to the time scale determined by the characteristic speeds of the system. In this second case the 
hyperbolic system with relaxation is said to be stiff and typically its solutions are well approximated 
by solutions of a suitable reduced set of conservation laws called equilibrium system |13j . 

Usually it is extremely difficult, if not impossible, to split the problem in separate regimes and 
to use different solvers in the stiff and non stiff regions. Thus one has to use the original relaxation 
system in the whole computational domain. The construction of schemes that work for all ranges 
of the relaxation time, using coarse grids that do not resolve the small relaxation time, has been 
studied mainly in the context of upwind methods using a method of lines approach combined with 
suitable operator splitting techniques [10, 25] and more recently in the context of central schemes 
[301135]. 

Splitting methods have been widely used for such problems. They are attractive because of 
their simplicity and robustness. Strang splitting provides second order accuracy if each step is at 
least second order accurate [55] • This property is maintained under fairly mild assumptions even 
for stiff problems |24) . However, Strang splitting applied to hyperbolic systems with relaxation 
reduces to first order accuracy when the problem becomes stiff. The reason is that the kernel of 
the relaxation operator is non trivial, which corresponds to a singular matrix in the linear case, 
and therefore the assumptions in [24] are not satisfied. 

Furthermore with a splitting strategy it is difficult to obtain higher order accuracy even in non 
stiff regimes (high order splitting schemes can be constructed, see [T5], but they are seldom used 
because of stability problems) . 

Recently developed Runge-Kutta schemes overcome this difficulties, providing basically the 
same advantages of the splitting schemes, without the drawback of the order restriction [TU1 121)1 14"5] . 

In this paper we will present some recent results on the development of high order, underresolved 
Runge-Kutta time discretization for such systems. In particular, using the formalism of implicit- 
explicit (IMEX) Runge-Kutta schemes [H [SJ [3H1 QT] we derive new IMEX schemes up to order 
3 that are strong-stability-preserving (SSP) for the limiting system of conservation laws (such 
methods were originally referred to as total variation diminishing (TVD) methods, see [18 ] fT9 ] 144 ] ) . 

To this aim, we derive general conditions that guarantee the asymptotic preserving property, 
i.e. the consistency of the scheme with the equilibrium system and the asymptotic accuracy, i.e. 
the order of accuracy is maintained in the stiff limit. For a stability analysis of IMEX schemes we 
refer to [38] . 

The rest of the paper is organized as follows. In Section 2 we introduce the general structure 
of IMEX Runge-Kutta schemes. Section 3 is devoted to IMEX Runge-Kutta schemes applied to 
hyperbolic systems with relaxation. In Section 4 we describe space discretization obtained by 
conservative finite-volume and finite difference schemes. In Section 5 we present some numerical 
results concerning the accuracy of IMEX schemes when applied to a prototype hyperbolic system 
with relaxation. Finally in Section 6 we investigate the performance of the schemes in several 
realistic applications to shallow waters, traffic flows and granular gases. 

2 IMEX Runge-Kutta schemes 

An IMEX Runge-Kutta scheme consists of applying an implicit discretization to the source terms 
and an explicit one to the non stiff term. When applied to system ([T]) it takes the form 

[/« = U n - AtJ2^jd x F(U u) ) + AtJ^^j-RiU^), 
U n+1 = U n - At^2wid x F(U^) + At^2wi-R(U^). 

i=l i=l E 

The matrices A = (5^ ) , = for j > i and A — (ay ) are v x v matrices such that the resulting 
scheme is explicit in F, and implicit in R. An IMEX Runge-Kutta scheme is characterized by 
these two matrices and the coefficient vectors w — (u>i, . . . , w 1/ ) T , w — (wi, . . . , w 1/ ) T . 

Since the simplicity and efficiency of solving the algebraic equations corresponding to the im- 
plicit part of the discretization at each step is of paramount importance, it is natural to consider 

2 



(3) 
(4) 



diagonally implicit Runge-Kutta (DIRK) schemes 22 for the source terms (a^ = 0, for j > i). 

IMEX Runge-Kutta schemes can be represented by a double tableau in the usual Butcher 
notation, 



w 



T 

W 



where the coefficients c and c used for the treatment of non autonomous systems, are given by the 
usual relation 



i-i 



Cj — ^ fl'ij, Cj — ^ aij. (5) 
i=i j'=i 

The use of a DIRK scheme for i? is a sufficient condition to guarantee that F is always evaluated 
explicitly. 

In the case of system ©, in order to obtain a numerical scheme, a suitable space discretization 
of equations ©-(01 is required. This discretization can be performed at the level of the original 
system © in which case one has a system of ODEs that is then solved in time by the IMEX scheme 
(method of lines). Alternatively one can apply a suitable space discretization directly to the time 
discrete equations ©-©. 

Finally we remark that previously developed schemes, such as the Additive semi-implicit Runge- 
Kutta methods of Zhong [49], and the splitting Runge-Kutta methods of Jin et al. [26], [10] can 
be rewritten as IMEX-RK schemes l38l. 



2.1 Order conditions 

The general technique to derive order conditions for Runge-Kutta schemes is based on the Taylor 
expansion of the exact and numerical solution. 

In particular, conditions for schemes of order p are obtained by imposing that the solution of 
system ([5]) at time t = to + Ai, with a given initial condition at time to, agrees with the numerical 
solution obtained by one step of a Runge-Kutta scheme with the same initial condition, up to order 
At p+1 . 

Here we report the order conditions for IMEX Runge-Kutta schemes up to order p — 3, which 
is already considered high order for PDEs problems. 

We apply scheme ©-© to system ©, with e = 1. We assume that the coefficients Cj, Cj, a^, 
aij satisfy conditions ©. Then the order conditions are the following 

First order 

^u>i = l, ^2wi = l. (6) 

8=1 1=1 

u>iCi = 1/2, ^ Wid = 1/2, (7) 

i i 

J2 ^ = 1/2, J2 w ^ = l i 2 - ( 8 ) 

i i 

Third order 

y^ WjaijCj = 1/6, ^2 WiCiCi = 1/3, y]wjaijCj = 1/6, ^WiC;C; = 1/3, (9) 

ij i ij i 

J2 i3 WiaijCj = 1/6, V, , <r,a,,c, = 1/6, V, ; ir,a,,c, = 1/6, 
J2ij WiCiijCj = 1/6, J^ij w i<ii]C] = !/6, w iai]Cj = 1/6, (10) 

J2i mciCi = 1/3, J2i WiCiCi = 1/3, J^i WiCiCi = 1/3, £\ lyjCjCj = 1/3. 



Second order 
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IMEX-RK 


Number of coupling conditions 


order 


General case 


Wi = Wi 


c = c 


c = c and it)j = u>; 


1 














2 


2 











3 


12 


3 


2 





4 


56 


21 


12 


2 


5 


252 


110 


54 


15 


6 


1128 


528 


218 


78 



Table 1: Number of coupling conditions in IMEX Runge-Kutta schemes 

Conditions ©, (J7J, © are the standard order conditions for the two tableau, each of them taken 
separately. Conditions ([SJ and (fTOf are new conditions that arise because of the coupling of the 
two schemes. 

The order conditions will simplify a lot if c = c. For this reason only such schemes are considered 
in [3]. In particular, we observe that, if the two tableau differ only for the value of the matrices A, 
i.e. if Cj = c, and Wi = Wi, then the standard order conditions for the two schemes are enough to 
ensure that the combined scheme is third order. Note, however, that this is true only for schemes 
up to third order. 

Higher order conditions can be derived as well using a generalization of Butcher 1-trees to 
2-trees, see |11) . However the number of coupling conditions increase dramatically with the order 
of the schemes. The relation between coupling conditions and accuracy of the schemes is reported 
in Table [U 



3 Applications to hyperbolic systems with relaxation 

In this section we give sufficient conditions for asymptotic preserving and asymptotic accuracy 
properties of IMEX schemes. This properties are strongly related to L-stability of the implicit part 
of the scheme. 

3.1 Zero relaxation limit 

Let us consider here one-dimensional hyperbolic systems with relaxation of the form ©. The 
operator R : M. N — s> R N is called a relaxation operator, and consequently © defines a relaxation 
system, if there exists a constant n x N matrix Q with rank(Q) = n < N such that 

QR(U) = \/UeR N . (11) 

This gives n independent conserved quantities u = QU. Moreover we assume that equation 
R(U) — can be uniquely solved in terms of u, i.e. 

U = £{u) such that R(£(u)) = 0. (12) 

The image of £ represents the manifold of local equilibria of the relaxation operator R. 

Using (jlip in ([2]) we obtain a system of n conservation laws which is satisfied by every solution 
of © 

d t (QU)+d x (QF(U)) = 0. (13) 

For vanishingly small values of the relaxation parameter e from © we get R(U) — which by 
(fT"2"|) implies U = £(u). In this case system © is well approximated by the equilibrium system 

d t u + d x G{u) = 0, (14) 

where G(u) = QF{£{u)). 

System (fl4"|) is the formal limit of system ([lj as e — s> 0. The solution u(x, t) of the system will be 
the limit of QU, with U solution of system ([1]), provided a suitable condition on the characteristic 
velocities of systems (HJ and (fT4"]) is satisfied (the so called subcharacteristic condition, see [4"8"lll3j .) 
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3.2 Asymptotic properties of IMEX schemes 

We start with the following 

Definition 3.1 We say that an IMEX scheme for system (0) in the form fJp-Q) is asymptotic 
preserving (AP) if in the limit e — > the scheme becomes a consistent discretization of the limit 
equation 

Note that this definition does not imply that the scheme preserves the order of accuracy in t in 
the stiff limit e — > 0. In the latter case the scheme is said asymptotically accurate. 

In order to give sufficient conditions for the AP and asymptotically accurate property, we make 
use of the following simple 

Lemma 3.1 If all diagonal elements of the triangular coefficient matrix A that characterize the 
DIRK scheme are non zero, then 

MmR(U {l) )={). (15) 

e->0 

Proof: 

In the limit e — > from ([3]) we have 

i 

^ aij R{W) = Q, t = l,...,i/. 
Since the matrix A is non-singular, this implies R(U l ) = 0, i = 1, . . . , v. 

■ 

In order to apply the previous lemma, the vectors of c and c cannot be equal. In fact c\ = 
whereas c\ ^ 0. Note that if c\ = but a« ^ for i > 1, then we still have lim e _j.o R(U' % ') = 
for i > 1 but lim e _j.o R(U^) ^ in general. The corresponding scheme may be inaccurate if the 
initial condition is not "well prepared" (R(Uq) / 0). In this case the scheme is not able to treat 
the so called "initial layer" problem and degradation of accuracy in the stiff limit is expected (see 
Section 5 and references [1011361135].) On the other hand, if the initial condition is "well prepared" 
(i?(C/ (0) ) = 0), then relation ([15]), i = l,...,v holds even if an = c\ = 0. 

Next we can state the following 

Theorem 3.1 If det A ^ in the limit e — » 0, the IMEX scheme (0)-|^]) applied to system (0) 
becomes the explicit RK scheme characterized by {A, w, c) applied to the limit equation |J^[ ). 

Proof: 

As in the case of the continuous system, multiplying equations ([3]) and ([4]) by the matrix Q, and 
making use of the relation QR(U) = 0VJ7, we obtain 

i-i 

u« =uo + h^OijQFiU®), (16) 

and 

Ul =u + hJ2wiQF(UW), (17) 
i=l 

where 

U W = QC7«, i = l,...v, u = QU , ui = QU x . 
From the previous lemma, in the stiff limit one has 

R(U (i) ) =0,i = l,...v, 

and, from property (|12p of the relaxation operator R, this implies 

U {l) =£(u^), i = l,...i/. 
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Substituting this expression in (|16l) - p7|) one obtains 

i-l 

U W =u + hJ2aijG(u^), (18) 

3=1 

and 

V 

u 1 =u Q + h y £ j w i G{u®), (19) 

i=l 

where = QF(S(u)). 

■ 

Clearly one may claim that if the implicit part of the IMEX scheme is A-stable or L-stablc 
the previous theorem is satisfied. Note however that this is true only if the tableau of the implicit 
integrator does not contain any column of zeros that makes it reducible to a simpler A-stable or 
L-stable form. Some remarks are in order. 

Remarks: 

i) There is a close analogy between hyperbolic systems with stiff relaxation and differential 
algebraic equations (DAE) [3]. The limit system as e — »■ is the analog of an index 1 DAE, 
in which the algebraic equation is explicitly solved in terms of the differential variable. In the 
context of DAE, the initial condition that we called "well prepared" is called "consistent" . 

ii) This result does not guarantee the accuracy of the solution for the N — n non conserved 
quantities. In fact, since the very last step in the scheme is not a projection towards the local 
equilibrium, a final layer effect occurs. The use of stiffly accurate schemes (i.e. schemes for 
which a u j = Wj,j = 1, . . . , ii) in the implicit step may serve as a remedy to this problem. 

iii) The theorem guarantees that in the stiff limit the numerical scheme becomes the explicit RK 
scheme applied to the equilibrium system, and therefore the order of accuracy of the limiting 
scheme is greater or equal to the order of accuracy of the original IMEX scheme. 

When constructing numerical schemes for conservation laws, one has to take a great care in 
order to avoid spurious numerical oscillations arising near discontinuities of the solution. This is 
avoided by a suitable choice of space discretization (see Section 4) and time discretization. 

Solution of scalar conservation equations, and equations with a dissipative source have some 
norm that decreases in time. It would be desirable that such property is maintained at a discrete 
level by the numerical method. If U n represents a vector of solution values (for example obtained 
from a method of lines approach in solving (fT4|) ) we recall the following [44] 

Definition 3.2 A sequence {U n } is said to be strongly stable in a given norm \ \ ■ \ \ provided that 
\\U n+1 \\ < \\U n \\ for all n > 0. 

The most commonly used norms are the Tl/-norm and the infinity norm. A numerical scheme 
that maintains strong stability at discrete level is called Strong Stability Preserving (SSP). 

Here we review some basic facts about RK-SSP schemes. First, it has been shown [19] under 
fairly general conditions that high order SSP schemes are necessarily explicit. Second, observe that 
a generic explicit RK scheme can be written as 

1/(0) = V n 
i-l 

= Y.^ u(k) +&tp ik L(uW), i = l,...,u, (20) 

fc=0 

U n + l = u( v ) 

where onk > and otik = only if = 0. This representation of RK schemes (which is not 
unique) can be converted to a standard Butcher form in a straightforward manner. Observe that 
for consistency, one has 5Z1 = q ctik = 1- It follows that if the scheme can be written in the form 



G 



(120[) with non negative coefficients ft/- then it is a convex combination of Forward Euler steps, with 
step sizes ftfe/a^At. A consequence of this is that if Forward Euler is SSP for At < At* , then the 
RK scheme is also SSP for At < cAt* , with c = mm lk (a ik / /3 ik ) [iUllTO] . 

The constant c is a measure of the efficiency of the SSP-RK scheme, therefore for the appli- 
cations it is important to have c as large as possible. For a detailed description of optimal SSP 
schemes and their properties see |44j . 

By point iii) of the above remarks it follows that if the explicit part of the IMEX scheme is 
SSP then, in the stiff limit, we will obtain an SSP method for the limiting conservation law. This 
property is essential to avoid spurious oscillations in the limit scheme for equation p4[) . 

In this section we give the Butcher tableau for the new second and third order IMEX schemes 
that satisfy the conditions of Theorem 13.11 In all these schemes the implicit tableau corresponds 
to an L-stable scheme, that is w T A~ 1 e — 1, e being a vector whose components are all equal to 1 
[22] . whereas the explicit tableau is SSP/c, where k denotes the order of the SSP scheme. Several 
examples of asymptotically SSP schemes are reported in tables We shall use the notation 

SSP£;(s, cr,p), where the triplet (s, <7,p) characterizes the number s of stages of the implicit scheme, 
the number a of stages of the explicit scheme and the order p of the IMEX scheme. 












7 


7 





1 


1 





1-7 


I-27 


7 




1/2 


1/2 
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1/2 



Table 2: Tableau for the explicit (left) implicit (right) IMEX-SSP2(2,2,2) L-stable scheme 
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Table 3: Tableau for the explicit (left) implicit (right) IMEX-SSP2 (3,2,2) stiffly accurate scheme 
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1/3 


1/3 


1/3 




1/3 


1/3 


1/3 




1/3 


1/3 


1/3 



Table 4: Tableau for the explicit (left) implicit (right) IMEX-SSP2(3,3,2) stiffly accurate scheme 

4 IMEX-WENO schemes 

For simplicity we consider the case of the single scalar equation 

u t + f(u) x = -r(u). (21) 

We have to distinguish between schemes based on cell averages (finite volume approach, widely 
used for conservation laws) and schemes based on point values (finite difference approach). 
Let Ax and At be the mesh widths. We introduce the grid points 

xj=jAx, Xj+i/2 = Xj + -Ax, j = ...,-2,-1,0,1,2,... 
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2/3 



Table 5: Tableau for the explicit (left) implicit (right) IMEX-SSP3(3,3,2) L-stable scheme 
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1/6 


1/6 


2/3 







1/6 


1/6 


2/3 



a = 0.24169426078821, /3 = 0.06042356519705 r/ = 0.12915286960590 
Table 6: Tableau for the explicit (left) implicit (right) IMEX-SSP3(4,3,3) L-stable scheme 

and use the standard notations 

i r x i+i/2 

itf = u(xj,t n ), u? = — / u(x,t n )dx. 
4.1 Finite volumes 

Integrating equation (|2ip on Ij = [ac^-i/a, £3+1/2] and dividing by Aa; we obtain 

1ft = -^[/M%+V2>*)) - /(«(*i-l/2,«))] + ^H,- (22) 

In order to convert this expression into a numerical scheme, one has to approximate the right hand 
side with a function of the cell averages {u(t)}j, which are the basic unknowns of the problem. 

The first step is to perform a reconstruction of the unknown function u{x, t) by a piecewise 
polynomial, starting from cell averages Uj(t). Given {uj}, compute a piecewise polynomial recon- 
struction 

3 

where Pj(x) is a polynomial satisfying some accuracy and non oscillatory property, and Xj( x ) is 
the indicator function of cell Ij. For first order schemes, the reconstruction is piecewise constant, 
while second order schemes can be obtained by a piecewise linear reconstruction. Higher order 
schemes are obtained by higher order polynomials. 

The reconstruction step is very delicate, because the function u(x, t) may be discontinuous. To 
this goal, suitable techniques, such as Weighted Essentially Non Oscillatory (WENO), have been 
developed. The reader can consult the review by Shu [4TJ and references therein for a more detailed 
description of high order non oscillatory reconstructions. 

The flux function at the edge of the cell can be computed by using a suitable numerical flux 
function, consistent with the analytical flux, 

f(u(x j+1/2 ,t)) F(uJ +1/2 ,u+ +1/2 ), 

where the quantities u^,, ,„ are obtained from the reconstruction, as u^,, , = lim_._± IZ(x). 
Example of numerical flux functions are the Godunov flux 

F( U j+l/2> U j+l/2) = /( U ( U j+l/2> U j+l/2))> 
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where u *( u J+1 / 2 ' ^+1/2) * s ^ ne som tion of the Riemann problem at Xj +1 / 2 , corresponding to the 



states u • 1 1/2 and u + , , and the Local Lax Friedrichs flux (also known as Rusanov flux), 



F ( u j+i/2> 4+1/2) = \\U( u 3 +i/2) + /(«i +1 /a)) ~ Q <ii/2 - Vi/a'l' 

where a = max„, |/'(w)|, and the maximum is taken over the relevant range of w. 

The two examples constitute two extreme cases of numerical fluxes: the Godunov flux is the 
most accurate and the one that produces the best results for a given grid size, but it is very 
expensive, since it requires the solution to the Riemann problem. Local Lax-Friedrichs flux, on the 
other hand, is less accurate, but much cheaper. This latter is the numerical flux that has been used 
throughout all the calculations performed for this paper. The difference in resolution provided by 
the various numerical fluxes becomes less relevant with the increase in the order of accuracy of the 
method. 

The right hand side of Eq. (f22|) contains the average of the source term r (u) instead of the source 
term evaluated at the average of u, r(u). The two quantities agree within second order accuracy 



r(u) j = r(%) + 0(Ax 2 ). 

This approximation can be used to construct schemes up to second order. 

First order (in space) semidiscrete schemes can be obtained using the numerical flux function 
F(uj,u j+ i) in place of f(u(x j+1 / 2 , t)), 

duj _ F(uj,u j+1 ) - F{u _ u u 3 ) | 1 
dt Ax e 3 

Second order schemes are obtained by using a piecewise linear reconstruction in each cell, and 
evaluating the numerical flux on the two sides of the interface 

du F ( U 3 +l/2> u t+i/2^ ~ F ( u j-i/2> "i-i/a) 1 ,_ s 

—— = — ' — ' — — ' — ' h -r(ui) 

dt Ax e K °' 

The quantities at cell edges are computed by piecewise linear reconstruction. For example, 

Ax , 

U j+l/2 — U 3 + ~2~ U 3 

where the slope u'j is a first order approximation of the space derivative of u(x,t), and can be 
computed by suitable slope limiters (see, for example, |29j for a discussion on TVD slope limiters.) 

For schemes of order higher than second, a suitable quadrature formula is required to approxi- 
mate g(u)j. For example, for third and fourth order schemes, one can use Simpson's rule 



r(u) 3 f» -(r(uJ_ 1/2 ) + Ar{u 3 ) + r(u j+1/2 )) 

where the pointwise values uf^^, Uj, uj +1 ^ 2 are obtained from the reconstruction. 

For a general problem, this has the effect that the source term couples the cell averages of 
different cells, thus making almost impractical the use of finite volume methods for high order 
schemes applied to stiff sources. 

Note, however, that in many relevant cases of hyperbolic systems with relaxation the implicit 
step, thanks to the conservation properties of the system, can be explicitly solved, and finite volume 
methods can be successfully used. We mention here all relaxation approximation of Jin-Xin type 
|28j . some simple discrete velocity models, such as Carlemann and Broadwell models, monoatomic 
gas in Extended Thermodynamics, semiconductor models, and shallow water equations. 

4.2 Finite differences 

In a finite difference scheme the basic unknown is the pointwise value of the function, rather than 
its cell average. Osher and Shu observed that it is possible to write a finite difference scheme in 
conservative form |42j . Let us consider the equation 

du df 1 . 
dt + d~x = S r{u) ' 
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and write 



df ( , + ,, f(u(x+^,t))-f(u(x-^,t)) 



dx ' Ax 
The relation between / and / is the following. Let us consider the sliding average operator 



1 f * 

u(x,t) =< u > x = — / u(£,t)d£. 
Ax J x _^l 



Differentiating with respect to x one has 

du 1 . . Ax . . Ax 
Tx = A^ {U{X + — >*)-<*- ~2 '*))■ 

Therefore the relation between / and / is the same that exists between u(x, t) and u(x, t), namely, 
ffux function / is the cell average of the function /. This also suggests a way to compute the flux 
function. The technique that is used to compute pointwise values of u(x,t) at the edge of the cell 
from cell averages of u can be used to compute f{u(xj + i/ 2 , t)) from f(u(xj, t)). This means that in 
the finite difference method it is the flux function which is computed at Xj and then reconstructed 
at Xj+i/2- But the reconstruction at 2j+i/2 may be discontinuous. Which value should one use? 
A general answer to this question can be given if one considers flux functions that can be split 

f(u) = f+(u) + f-(u), (24) 

with the condition that 

df+(u) df-(u) , s 

J , >0, J —r < 0- 25 
du du 

There is a close analogy between flux splitting and numerical flux functions. In fact, if a flux can 
be split as ([2"4f. then 

F(a,b) = f+(a) + f-(b) 



will define a monotone consistent flux, provided condition (|25[) is satisfied. Together with non 
oscillatory reconstructions (such as WENO) and SSP time discretization, the monotonicity condi- 
tion will ensure that the overall scheme will not produce spurious numerical oscillations (see, for 
example, [29| [41].) 

This is the case, for example, of the local Lax-Friedrichs flux. A finite difference scheme 
therefore takes the following form 



dui 1 . » - 1 



F j+i/2 = f + (u j+1/2 ) + f («J +1/3 ); 

f + ( u J + i/ 2 ) i s obtained by 

• computing f + (ui) and interpret it as cell average of / + ; 

• performing pointwise reconstruction of / + in cell j, and evaluate it in Xj+i/%. 
/'("jlirt) i s obtained by 

• computing f~(ui), interpret as cell average of /~; 

• performing pointwise reconstruction of /~ in cell j + 1, and evaluate it in lEj+i/a- 

A detailed account on high order finite difference schemes can be found in [41) . 

At variance with finite volume schemes, since the source is evaluated pointwise, finite difference 
schemes do not couple the cells. This is a big advantage in the cases in which the implicit step can 
not be explicitly solved. 
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Remarks: 



i) Finite difference can be used only with uniform (or smoothly varying) mesh. In this respect 
finite volume are more flexible, since they can be used even with unstructured grids. 

ii) The treatment presented for the scalar equation can be extended to systems, with a minor 
change of notation. In particular, the parameter a appearing in the local Lax-Friedrichs flux 
will be computed using the spectral radius of the Jacobian matrix. 

iii) For schemes applied to systems, better results are usually obtained if one uses characteristic 
variables rather than conservative variables in the reconstruction step. The use of conser- 
vative variables may result in the appearance of small spurious oscillations in the numerical 
solution. For a treatment of this effect see, for example, [39] . 



5 Numerical tests 

In this section we investigate numerically the convergence rate and the zero relaxation limit behav- 
ior of the schemes. To this aim we apply the IMEX-WENO schemes to the Broadwell equations 
of rarefied gas dynamics [10, 26, 30, 35 . In all the computations presented in this paper we used 
finite difference WENO schemes with Lax-Friedrichs flux and conservative variables. Of course the 
sharpness of the resolution of the numerical results can be improved using a less dissipative flux. 

As a comparison, together with the new IMEX-SSP schemes, we have considered the second 
order ARS(2,2,2) method presented in [3], which was the one recommended by the authors. For 
this scheme, since c\ — 0, there will be a degradation of accuracy due to the initial layer. 

In all figures, the value of N represents the number of grid points in space. 

The kinetic model is characterized by a hyperbolic system with relaxation of the form with 

U = (p, to, z), F(U) = (to, z, to), R(U) = (o, 0, ^{p 2 + to 2 - 2pz) 

Here e represents the mean free path of particles. The only conserved quantities are the density p 
and the momentum to. 

In the fluid-dynamic limit e — > we have 

p 2 + to 2 , . 

z = z E ^ t— , 26 

2p 

and the Broadwell system is well approximated by the reduced system (|14[) with 

(I \ f j i 

which represents the corresponding Euler equations of fluid-dynamics. 
Convergence rates 

We have considered a periodic smooth solution with initial data as in [101 [30j given by 

t n \ * ■ 27ra , ^ 1 • 27rx / ^ p(x,0) 2 +m(x, 0) 2 . 
p(x,0) = l + a p sm—, v {x,0 = - +a v sm— , z x,0 ) = a / K '-> t ■ 27 
L 2 L 2p(x,0) 

In our computations we used the parameters 

a p = 0.3, a v = 0.1, a z — 1.0 (no initial layer) and a z = 0.2 (initial layer), L = 20, (28) 

and we integrate the equations for t 6 [0,30]. A Courant number At/ Ax — 0.6 has been used. 
The plots of the relative error are given in Figure [T] 

Notice how, in absence of initial an layer, all schemes tested have the prescribed order of 
accuracy both in the non stiff and in the stiff limit, with some degradation of the accuracy at 
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intermediate regimes. As expected, scheme ARS(2,2,2), shows a degradation of the accuracy when 
an initial layer is present. 

Next we test the shock capturing properties of the schemes in the case of non smooth solutions 
characterized by the following two Riemann problems [10] 

Pi = 2, mi = 1, zi = l, x < 0.2, ^ 

p r = 1, m r = 0.13962, z r = 1, x > 0.2, 

Pi = 1, m = 0, zi = 1, x < 0, 

= 0.2, TO,. =0, 2 r = 1, X > 0. 

For brevity we report the numerical results obtained with the second order IMEX-SSP2(2,2,2) and 
third order IMEX-SSP3(4,3,3) schemes that we will refer to as IMEX-SSP2-WENO and IMEX- 
SSP3-WENO respectively. The result are shown in Figures and |3] for a Courant number At/ Ax = 
0.5. Both schemes, as expected, give an accurate description of the solution in all different regimes 
also using coarse meshes that do not resolve the small scales. In particular the shock formation in 
the fluid limit is well captured without spurious oscillations. We refer to [TU1 HE1 [301 ISS1 12] for a 
comparison of the present results with previous ones. 

6 Applications 

Finally we present some numerical results obtained with IMEX-SSP2-WENO and IMEX-SSP3- 
WENO concerning situations in which hyperbolic systems with relaxation play a major role in 
applications. The results have been obtained with N = 200 grid points. As usual the reference 
solution is computed on a much finer grid. 

6.1 Shallow water 

First we consider a simple model of shallow water flow [26] 

d t h + d x (hv) = 0, 

d t (hv)+d x {h+\h 2 ) = \(\-v)> 

where h is the water height with respect to the bottom and hv the flux. 

The zero relaxation limit of this model is given by the inviscid Burgers equation. 
The initial data we have considered is |26j 

h 2 

h = 1 + 0.2sin(87re), hv = —, (32) 

with x <E [0,1]. The solution at t = 0.5 in the stiff regime e = 10~ 8 using periodic boundary 
conditions is given in Figure [4] For IMEX-SSP2-WENO the dissipative effect due to the use of 
the Lax-Friedrichs flux is very pronounced. As expected this effect becomes less relevant with the 
increase of the order of accuracy. We refer to [261 for a comparison with the present results. 

6.2 Traffic flows 

In [5] a new macroscopic model of vehicular traffic has been presented. The model consists of a 
continuity equation for the density p of vehicles together with an additional velocity equation that 
describes the mass flux variations due to the road conditions in front of the driver. The model can 
be written in conservative form as follows 

d t p + d x (pv) = 0, 
d t (pw) + d x (vpw) = A^(V(p)-v), 
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where w — v + P(p) with P{p) a given function describing the anticipation of road conditions in 
front of the drivers and V(p) describes the dependence of the velocity with respect to the density 
for an equilibrium situation. The parameter T is the relaxation time and A > is a positive 
constant. 

If the relaxation time goes to zero, under the subcharacteristic condition 

-P'{p) < V'{p) < 0, P >0, 
we obtain the Lighthill-Whitham [?5] model 

d t p + d x (pV(p))=0, (34) 
A typical choice for the function P(p) is given by 

P(p) = 

where p m is a given maximal density and c v a constant with dimension of velocity. 

In our numerical results we assume A — 1 and an equilibrium velocity V(p) fitting to experi- 
mental data [7] 

V2 + arctan(q£fe^) 

Vi -P>- Vm V2 + arctan(a/3) 

with a = 11, /3 = 0.22 and v m a maximal speed. We consider 7 = and, in order to fulfill the 
subcharacteristic condition, assume c v — 2. All quantities are normalized so that v m — 1 and 
p m = 1- 

We consider a Riemann problem centered at x — with left and right states 

PL = 0.05, v L = 0.05, p R = 0.05, v R = 0.5. (35) 

The solution at t = 1 for T = 0.2 is given in Figure [5] The figure shows the development of the 
density of the vehicles. Both schemes gives very similar results. Again, in the second order scheme 
the shock is smeared out if compared to the third order case. See [7] for more numerical results. 

6.3 Granular gases 

We consider the continuum equations of Euler type for a granular gas [2S1 27] . These equations 
have ben derived for a dense gas composed of inelastic hard spheres. The model reads 

Pt + (pu)x = 0, 
(pu) t + {pu 2 + p) x = pg, ^ 

(lpu 2 + lpT^+(\pu* + lupT + P u^ = -ii-£!l G (pV 2 T 3 / 2 , 

where e is the coefficient of restitution, g the acceleration due to gravity, e a relaxation time, p is 
the pressure given by 

p = pT(l + 2(l + e)G(p)), 
and G(p) is the statistical correlation function. In our experiments we assume 

G( P ) = Jl-(— 

where v = <j 3 pir/6 is the volume fraction, a is the diameter of a particle, and vm — 0.64994 is 3D 
random close-packed constant. 

We consider the following initial data [32] on the interval [0, 10] 

p = 34.37746770, v = 18, P = 1589.2685472, (37) 
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which corresponds to a supersonic flow at Mach number M a — 7 (the ratio of the mean fluid 
speed to the speed of sound). Zero- flux boundary condition have been used on the bottom (right) 
boundary whereas on the top (left) we have an ingoing flow characterized by (|37[) . 

The values of the restitution coefficient and the particle diameter have been taken e = 0.97 
and a — 0.1. We report the solution at t = 0.2 with e = 0.01 in Figure [6] (see [32] for similar 
results). Due to the nonlinearity of the source term the implicit solver has been efficiently solved 
using Newton's method in each cell thanks to the use of finite difference space discretization. 

Both methods provide a good description of the shock that propagates backward after the 
particles impact with the bottom. Note that the second order method provides excessive smearing 
of the layer at the right boundary. This problem is not present in the third order scheme. However 
due to the use of conservative variables we can observe the presence of small spurious oscillations 
in the pressure profile. 
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Relative error in density vs N, e = 1 , no initial layer 



Relative error in density vs N, e = 1 , with initial layer 
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Relative error in density vs N, e = 10 , no initial Is 



Relative error in density vs N, e = 10 , with initial layer 
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Figure 1: Relative errors for density p in the Broadwell equations with initial data ([28]) . Left 
column a z — 1.0 (no initial layer), right column a z = 0.2 (initial layer). From top to bottom, 
e = 1.0, 10" 3 , 10" 6 . 
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Figure 2: Numerical solution of the Broadwell equations with initial data (150)) for p(o), m(*) and 
z(+) at time i = 0.5. Left column IMEX-SSP2-WENO scheme, right column IMEX-SSP3-WENO 
scheme. From top to bottom, e = 1.0, 0.02, 10~ 8 . 
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IMEX-SSP2-WEN0, e-10 , 1=0.25, N=200 



IMEX-SSP3-WENO,e=10 , 1=0.25, N=200 
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Figure 3: Numerical solution of the Broadwell equations with initial data (flTT]) for p(o), m(*) and 
z{+) at time t = 0.25 for e = 10~ 8 . Left IMEX-SSP2-WENO scheme, right IMEX-SSP3-WENO 
scheme. 
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Figure 4: Numerical solution of the shallow water model with initial data (|32| for h(o) and hv(*) 
at time t = 0.5 for e = 10" 8 . Left IMEX-SSP2-WENO scheme, right IMEX-SSP3-WENO scheme. 



IMEX-SSP2-WENO, t=0.2. t=l N=200 
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Figure 5: Numerical solution of the traffic model with initial data (|35|) for p(o) and pv(*) at time 
t = 1 for e = 0.2. Left IMEX-SSP2-WENO scheme, right IMEX-SSP3-WENO scheme. 
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IMEX-SSP2-WEN0. £=0.01 , t=0.2, N=200 



IMEX-SSP3-WENO, 6=0.01 , t=0.2, N=200 




Figure 6: Numerical solution of the hydrodynamical model of a granular gas with initial data (|37|) . 
Left column IMEX-SSP2-WENO scheme, right column IMEX-SSP3-WENO scheme. From top to 
bottom, mass fraction v, velocity pu and pressure p. 
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